%**************************************************************************************************
%Author: Yong Wang
%Last Edited: July 1, 2010
%Email: ywang@csu.edu.cn
%Reference: Differential Evolution with Composite Trial Vector Generation Strategies
%                                                    and Control Parameters
%                           IEEE Transactions on Evolutionary Computation, Accepted
%**************************************************************************************************
function [position,value,Pg]= CoDE(~,popsize,dimension,xmax,xmin,vmax,vmin,MaxDT,benchmark_func,func_num,~)

    max_FES=10000*dimension;

    problem = func_num;

    % Define the dimension of the problem
    n = dimension;
    lu = [xmin * ones(1, n); xmax * ones(1, n)];

    % Main body'

    fitold=1;

    Pg=zeros(1,max_FES);    %迭代过程中的最优适应度值 

    % Initialize the main population
    p = repmat(lu(1, :), popsize, 1) + rand(popsize, n) .* (repmat(lu(2, :) - lu(1, :), popsize, 1));

    % Evaluate the objective function values
    for i = 1:popsize
        fit = benchmark_func(p', problem);
    end
    fit=fit';
    % Record the number of function evaluations (FES)
    FES = popsize;
    m=0;
    while FES < max_FES
        m=m+1;
        pTemp = p;
        fitTemp = fit;
        
        % uSet: the set of trial vectors
        uSet = zeros(3 * popsize, n);

        for i = 1 : popsize

            % the three control parameter settings
            F    = [1.0 1.0 0.8];
            CR = [0.1 0.9 0.2];

            % Uniformly and randomly select one of the control
            % parameter settings for each trial vector generation strategy
            paraIndex = floor(rand(1, 3) * length(F)) + 1;

            % Generate the trail vectors
            u = generator(p, lu, i, F, CR, popsize, n, paraIndex);

            uSet(i * 3 - 2 : 3 * i, :) = u;

            FES = FES + 3;

        end

        % Evaluate the trial vectors
        fitSet = benchmark_func(uSet', problem);
        fitSet=fitSet';
        for i = 1 : popsize

            % Choose the best trial vector from the three trial vectors
            [~, minID] = min(fitSet(3 * i - 2 : 3 * i, :));
            bestInd = uSet(3 * (i - 1) + minID, :);
            bestIndFit = fitSet(3 * (i - 1) + minID, :);
            % Choose the better one between the trial vector and the
            % target vector
            if fit(i) >= bestIndFit
                pTemp(i, :) = bestInd;
                fitTemp(i, :) = bestIndFit;

            end
        end
        p = pTemp;
        fit = fitTemp;


        [value,id]=min(fit);
        position=p(id,:);
   
        Pg(1,fitold:FES)=value;

        if mod(FES,1e4)<41
            fprintf("CoDE 第%d次评价，最佳适应度 = %e\n",FES,value);
        end

        fitold=FES;
        if FES>=max_FES
           break;
        end
    end
    
 
function u = generator(p, lu, i, F, CR, popsize, n, paraIndex)

%.... "rand/1/bin" strategy ....%

% Choose the indices for mutation
indexSet = 1 : popsize;
indexSet(i) = [];

% Choose the first Index
temp = floor(rand * (popsize - 1)) + 1;
index(1) = indexSet(temp);
indexSet(temp) = [];

% Choose the second index
temp = floor(rand * (popsize - 2)) + 1;
index(2) = indexSet(temp);
indexSet(temp) = [];

% Choose the third index
temp = floor(rand * (popsize - 3)) + 1;
index(3) = indexSet(temp);

% Mutation
v1 = p(index(1), :) + F(paraIndex(1)) .* (p(index(2), :) - p(index(3), :));

% Handle the elements of the mutant vector which violate the boundary
vioLow = find(v1 < lu(1, :));
if ~isempty(vioLow)
    v1(1, vioLow) = 2 .* lu(1, vioLow) - v1(1, vioLow);
    vioLowUpper = find(v1(1, vioLow) > lu(2, vioLow));
    if ~isempty(vioLowUpper)
        v1(1, vioLow(vioLowUpper)) = lu(2, vioLow(vioLowUpper));
    end
end

vioUpper = find(v1 > lu(2, :));
if ~isempty(vioUpper)
    v1(1, vioUpper) = 2 .* lu(2, vioUpper) - v1(1, vioUpper);
    vioUpperLow = find(v1(1, vioUpper) < lu(1, vioUpper));
    if ~isempty(vioUpperLow)
        v1(1, vioUpper(vioUpperLow)) = lu(1, vioUpper(vioUpperLow));
    end
end

% Binomial crossover
j_rand = floor(rand * n) + 1;
t = rand(1, n) < CR(paraIndex(1));
t(1, j_rand) = 1;
t_ = 1 - t;
u(1, :) = t .* v1 + t_ .* p(i, :);

                                %... "current to rand/1" strategy ...%

% The mechanism to choose the indices for mutation is slightly different from that of the classic
% "current to rand/1", we found that using the following mechanism to choose the indices for
% mutation can improve the performance to certain degree
index(1: 3) = floor(rand(1, 3) * popsize) + 1;

% Mutation
v2 = p(i, :) + rand * (p(index(1), :) - p(i, :)) + F(paraIndex(2)) .* (p(index(2), :) - p(index(3), :));

% Handle the elements of the mutant vector which violate the boundary
vioLow = find(v2 < lu(1, :));
if ~isempty(vioLow)
    v2(1, vioLow) = 2 .* lu(1, vioLow) - v2(1, vioLow);
    vioLowUpper = find(v2(1, vioLow) > lu(2, vioLow));
    if ~isempty(vioLowUpper)
        v2(1, vioLow(vioLowUpper)) = lu(2, vioLow(vioLowUpper));
    end
end

vioUpper = find(v2 > lu(2, :));
if ~isempty(vioUpper)
    v2(1, vioUpper) = 2 .* lu(2, vioUpper) - v2(1, vioUpper);
    vioUpperLow = find(v2(1, vioUpper) < lu(1, vioUpper));
    if ~isempty(vioUpperLow)
        v2(1, vioUpper(vioUpperLow)) = lu(1, vioUpper(vioUpperLow));
    end
end

% Binomial crossover is not used to generate the trial vector under this
% condition
u(2, :) = v2;

%... "rand/2/bin" strategy ...%

% Choose the indices for mutation
indexSet = 1 : popsize;
indexSet(i) = [];

% Choose the first index
temp = floor(rand * (popsize - 1)) + 1;
index(1) = indexSet(temp);
indexSet(temp) = [];

% Choose the second index
temp = floor(rand * (popsize - 2)) + 1;
index(2) = indexSet(temp);
indexSet(temp) = [];

% Choose the third index
temp = floor(rand * (popsize - 3)) + 1;
index(3) = indexSet(temp);
indexSet(temp) = [];

% Choose the fourth index
temp = floor(rand * (popsize - 4)) + 1;
index(4) = indexSet(temp);
indexSet(temp) = [];

% Choose the fifth index
temp = floor(rand * (popsize - 5)) + 1;
index(5) = indexSet(temp);

% Mutation
% The first scaling factor is randomly chosen from 0 to 1
v3 = p(index(1), :) + rand .* (p(index(2), :) - p(index(3), :)) + F(paraIndex(3)) .* (p(index(4), :) - p(index(5), :));

% Handle the elements of the mutant vector which violate the boundary
vioLow = find(v3 < lu(1, :));
if ~isempty(vioLow)
    v3(1, vioLow) = 2 .* lu(1, vioLow) - v3(1, vioLow);
    vioLowUpper = find(v3(1, vioLow) > lu(2, vioLow));
    if ~isempty(vioLowUpper)
        v3(1, vioLow(vioLowUpper)) = lu(2, vioLow(vioLowUpper));
    end
end

vioUpper = find(v3 > lu(2, :));
if ~isempty(vioUpper)
    v3(1, vioUpper) = 2 .* lu(2, vioUpper) - v3(1, vioUpper);
    vioUpperLow = find(v3(1, vioUpper) < lu(1, vioUpper));
    if ~isempty(vioUpperLow)
        v3(1, vioUpper(vioUpperLow)) = lu(1, vioUpper(vioUpperLow));
    end
end

% Binomial crossover
j_rand = floor(rand * n) + 1;
t = rand(1, n) < CR(paraIndex(3));
t(1, j_rand) = 1;
t_ = 1 - t;
u(3, :) = t .* v3 + t_ .* p(i, :);
end
   
  

    if FES<max_FES
        Pg(FES+1:max_FES)=value;
    else
        if FES>max_FES
            Pg(max_FES+1:end)=[];
        end
    end 
   
end
